library(AIUQ)
Example 1: Simulate BM and get estimated parameters using BM
model
set.seed(1)
## Simulation
sim_bm = simulation()
show(sim_bm)
## Frame size: 200 200
## Number of time steps: 200
## Number of particles: 50
## Stochastic process: BM
## Variance of background noise: 20
## sigma_bm: 1
## Plot simulated particle trajectory
plot_traj(sim_bm)
## Plot intensity profile for different frames
plot_intensity(sim_bm@intensity, sz=sim_bm@sz) #first frame
plot_intensity(sim_bm@intensity, sz=sim_bm@sz,frame=10, color=T) #10th frame, color image
## AIUQ method: use BM as fitted model
sam = SAM(sim_object=sim_bm, model_name='BM')
show(sam)
## Fitted model: BM
## Number of q ring: 99
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
## True parameters in the model: 2
## Estimated parameters in the model: 2.010371
## True variance of background noise: 20
## Estimated variance of background noise: 19.90164
## Maximum log likelihood value: -18150726
## Akaike information criterion score: 36301654
## Plot shifted intensity in reciprocal space
plot_I_q = matrix(sam@I_q[,1], sam@sz[1],sam@sz[2]) #first frame
plot3D::image2D(abs(fftshift(plot_I_q)),main="Foruier tansformed intensity")




## User can select q range via AIUQ_thr or index_q
sam = SAM(sim_object=sim_bm, AIUQ_thr=c(0.99,0.6))
#Note: Default model_name is "BM", it's ok to not specify this argument if want to fit with BM model
show(sam)
## Fitted model: BM
## Number of q ring: 99
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## True parameters in the model: 2
## Estimated parameters in the model: 2.008843
## True variance of background noise: 20
## Estimated variance of background noise: 20.05226
## Maximum log likelihood value: -8065526
## Akaike information criterion score: 16131176
sam = SAM(sim_object=sim_bm, index_q_AIUQ=5:50)
show(sam)
## Fitted model: BM
## Number of q ring: 99
## Index of wave number selected: 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## True parameters in the model: 2
## Estimated parameters in the model: 2.015093
## True variance of background noise: 20
## Estimated variance of background noise: 20.02955
## Maximum log likelihood value: -6198715
## Akaike information criterion score: 12397526
Example 4: Load csv file with intensity profile being a matrix with
dim time by space*space (T_SS_mat) and est parameters using BM
model
## Load intensity profile from csv file
file_loc = "/Users/yue/Downloads/"
file_name = "intensity_record_BM_beta_0_02_B_40_len_100_frame_size_100"
intensity = data.table::fread(paste(file_loc,file_name,".csv",sep=""))
intensity = intensity[,-1]
## AIUQ method: use BM as fitted model
mindt = 1 #minimum lag time
pxsz = 1 #pixel size
sz=c(100,100) #frame size of the image
sam = SAM(intensity=intensity,mindt=mindt,pxsz=pxsz,sz=sz)
## start of another optimization, initial values: 0.3625398 3.5165
show(sam)
## Fitted model: BM
## Number of q ring: 49
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## True parameters in the model: NA
## Estimated parameters in the model: 0.02040957
## True variance of background noise: NA
## Estimated variance of background noise: 20.11201
## Maximum log likelihood value: -2102329
## Akaike information criterion score: 4204761
#sam@msd_est #estimated MSD
## DDM method: use fixed A and B
sam = SAM(intensity=intensity,mindt=mindt,pxsz=pxsz,sz=sz, method="DDM_fixedAB")
show(sam)
## Fitted model: BM
## Number of q ring: 49
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## True parameters in the model: NA
## Estimated parameters in the model: 8.336566
## True variance of background noise: NA
## Estimated variance of background noise: 20.35795
Example 5: Load tif file (SST_array) from real experiment and est
parameters using FBM model
## Load real experiment video
file_loc = "/Users/yue/Desktop/Research/2.SAM/AIUQ_package/data/"
file_name = "4pct_PVA_100nm_25C"
intensity = bioimagetools::readTIF(paste(file_loc,file_name,".tif",sep=""),as.is=TRUE)
## Plot intensity profile for different frames
intensity_str = 'SST_array' #intensity profile structure
plot_intensity(intensity, intensity_str=intensity_str) #first frame
plot_intensity(intensity, intensity_str=intensity_str, frame=15) #15th frame
## AIUQ method: fitting using FBM model
mindt = 0.0309 #minimum lag time
pxsz = 0.2930 #pixel size
model_name = "FBM"
sam = SAM(intensity=intensity,intensity_str=intensity_str,mindt=mindt, pxsz=pxsz, uncertainty=T, AIUQ_thr=c(0.99,0), model_name=model_name)
show(sam)
## Fitted model: FBM
## Number of q ring: 255
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
## True parameters in the model: NA
## Estimated parameters in the model: 0.634021 1.014347
## True variance of background noise: NA
## Estimated variance of background noise: 150.4655
## Maximum log likelihood value: -41045311
## Akaike information criterion score: 82090786
## Plot reference MSD vs estimated MSD with uncertainty
msd_truth = 0.7589*sam@d_input
plot_MSD(object=sam, msd_truth=msd_truth)



Example 6: Simulate FBM and compare estimated parameters using SAM
and DDM methods
set.seed(1)
## Simulation
sim_fbm = simulation(sz=100,len_t=100,model_name="FBM")
show(sim_fbm)
## Frame size: 100 100
## Number of time steps: 100
## Number of particles: 50
## Stochastic process: FBM
## Variance of background noise: 20
## (sigma_fbm, Hurst parameter): 2 0.3
## AIUQ method
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name)
show(sam)
## Fitted model: FBM
## Number of q ring: 49
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## True parameters in the model: 8 0.6
## Estimated parameters in the model: 8.02211 0.5969453
## True variance of background noise: 20
## Estimated variance of background noise: 15.05404
## Maximum log likelihood value: -2414150
## Akaike information criterion score: 4828404
# Plot true MSD vs estimated MSD
plot_MSD(object=sam,msd_truth=sam@msd_truth)
## DDM method with fixed A and B, using all q
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_fixedAB")
show(sam)
## Fitted model: FBM
## Number of q ring: 49
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## True parameters in the model: 8 0.6
## Estimated parameters in the model: 6.561253 0.8439145
## True variance of background noise: 20
## Estimated variance of background noise: 20.55181
## DDM method with fixed A and B, using user defined q range
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_fixedAB", index_q_DDM=3:40)
show(sam)
## Fitted model: FBM
## Number of q ring: 49
## Index of wave number selected: 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## True parameters in the model: 8 0.6
## Estimated parameters in the model: 7.04659 0.8025151
## True variance of background noise: 20
## Estimated variance of background noise: 20.55181
## Add line of estimated MSD from DDM method with fixed A and B and user defined q range
lines(log10(sam@d_input[-1]),log10(sam@msd_est[-1]),type='p',
col='red', pch=24, cex=0.8)
## DDM method with estimated A and B,using all q
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_estAB")
show(sam)
## Fitted model: FBM
## Number of q ring: 49
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## True parameters in the model: 8 0.6
## Estimated parameters in the model: 4.001553 0.9926345
## True variance of background noise: 20
## Estimated variance of background noise: 37.45107
## DDM method with estimated A and, using user defined q range
sam = SAM(sim_object=sim_fbm, model_name=sim_fbm@model_name, method="DDM_estAB", index_q_DDM=3:40)
show(sam)
## Fitted model: FBM
## Number of q ring: 49
## Index of wave number selected: 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## True parameters in the model: 8 0.6
## Estimated parameters in the model: 5.036552 0.9891863
## True variance of background noise: 20
## Estimated variance of background noise: 37.86241

Example 7: User defined MSD structure
set.seed(1)
## Simulation
sim_bm = simulation(sz=100,len_t=100,sigma_bm=0.5)
show(sim_bm)
## Frame size: 100 100
## Number of time steps: 100
## Number of particles: 50
## Stochastic process: BM
## Variance of background noise: 20
## sigma_bm: 0.5
## User defined MSD structure: function of parameters and
# vector of lag times
msd_fn = function(param, d_input){
MSD = param[1]*d_input+param[2]*d_input^2
}
# show MSD and MSD gradient with a simple example
theta = c(2,1)
d_input = 0:10
model_name = "user_defined"
MSD_list = get_MSD_with_grad(theta=theta,d_input=d_input,model_name=model_name,
msd_fn=msd_fn)
MSD_list$msd
## [1] 0 3 8 15 24 35 48 63 80 99 120
## AIUQ method: fitting using user_defined model
sam = SAM(sim_object=sim_bm, model_name=model_name, msd_fn=msd_fn, num_param=2)
show(sam)
## Fitted model: user_defined
## Number of q ring: 49
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## True parameters in the model: 0.5
## Estimated parameters in the model: 0.5047796 0.000272548
## True variance of background noise: 20
## Estimated variance of background noise: 20.08427
## Maximum log likelihood value: -2289361
## Akaike information criterion score: 4578827
Example 8: Extract empirical/modeled Dqt and ISF outside/within
SAM
set.seed(1)
## Simulation
sim_bm = simulation()
show(sim_bm)
## Frame size: 200 200
## Number of time steps: 200
## Number of particles: 50
## Stochastic process: BM
## Variance of background noise: 20
## sigma_bm: 1
## AIUQ method: use BM as fitted model
sam = SAM(sim_object=sim_bm, model_name='BM')
show(sam)
## Fitted model: BM
## Number of q ring: 99
## Index of wave number selected: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
## True parameters in the model: 2
## Estimated parameters in the model: 2.010371
## True variance of background noise: 20
## Estimated variance of background noise: 19.90164
## Maximum log likelihood value: -18150726
## Akaike information criterion score: 36301654
# Compute empirical/modeled Dqt and ISF after SAM class is constructed
par(mfrow = c(2, 2))
## Compute observed dynamic image structure function (Dqt)
Dqt = get_dqt(sam) #q index range is all qs without user specify
Dqt = get_dqt(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6, main="Observed dynamic image structure function")
## Compute empirical intermediate scattering function (ISF)
ISF = get_isf(sam) #q index range is all qs without user specify
ISF = get_isf(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Empirical intermediate scattering function")
## Compute modeled Dqt
modeled_Dqt = modeled_dqt(sam) #q index range is all qs without user specify
modeled_Dqt = modeled_dqt(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(modeled_Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6, main="Modeled dynamic image structure function")
## Compute modeled ISF
modeled_ISF = modeled_isf(sam) #q index range is all qs without user specify
modeled_ISF = modeled_isf(sam, index_q = 1:50)
matplot(sam@d_input[-1], t(modeled_ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Modeled intermediate scattering function")

# Compute empirical/modeled within SAM class
## AIUQ method: use BM as fitted model and compute empirical Dqt within class
sam = SAM(sim_object=sim_bm, model_name='BM', output_dqt=T)
matplot(sam@d_input[-1], t(sam@Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6,main="Observed dynamic image structure function")
## AIUQ method: use BM as fitted model and compute empirical Dqt,ISF within class
sam = SAM(sim_object=sim_bm, model_name='BM', output_isf=T) # Note: empirical ISF is calculated through observed Dqt. By setting output_isf is true, Dqt is also computed even without "output_dqt=T".
matplot(sam@d_input[-1], t(sam@Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6,main="Observed dynamic image structure function")
matplot(sam@d_input[-1], t(sam@ISF[-c(51:99),]), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Empirical intermediate scattering function")
matplot(sam@d_input[-1], t(sam@ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Empirical intermediate scattering function")

## AIUQ method: use BM as fitted model and compute ISF within class
sam = SAM(sim_object=sim_bm, model_name='BM', output_modeled_isf=T)
matplot(sam@d_input[-1], t(sam@modeled_ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6,main="Modeled intermediate scattering function")
## AIUQ method: use BM as fitted model and modeled Dqt,ISF within class
sam = SAM(sim_object=sim_bm, model_name='BM', output_modeled_dqt=T) # Note: Modeled Dqt is calculated through modeled ISF. By setting output_modeled_dqt is true, modeled_ISF is also computed even without "output_modeled_isf=T".
matplot(sam@d_input[-1], t(sam@modeled_ISF), xlab=expression(paste(Delta, "t",sep="")), ylab="ISF", type="l", lwd=0.6, main="Modeled intermediate scattering function")
matplot(sam@d_input[-1], t(sam@modeled_Dqt), xlab=expression(paste(Delta, "t",sep="")), ylab="Dqt", type="l", lwd=0.6,main="Modeled dynamic image structure function")
